!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!> \brief
!! Darcy flow in a fractured medium.
!!
!! \n
!!
!! This test case is described in <a
!! href="http://scitation.aip.org/getpdf/servlet/GetPDFServlet?filetype=pdf&id=SJOCE3000026000005001667000001&idtype=cvips&prog=normal">[Martin,
!! Jaffr&eacute;, Roberts, 2005]</a>. The boundary regions are
!! numbered as follows:
!! <ul>
!!   <li> bottom boundary, from left to right: indexes 1, 2, 3
!!   <li> right boundary: index 4
!!   <li> top boundary, from right to left: indexes 5, 6, 7
!!   <li> left boundary: index 8.
!! </ul>
!!
!! This module can be used in both when the fracture is represented
!! explicitly and when it is represented as a one-dimensional
!! manifold. The parameter <tt>fd</tt> represents the fracture
!! width, which must be always positive, while <tt>fdg</tt>
!! represents the <em>geometrical</em> fracture width, which can
!! be equal to <tt>fd</tt> (resolved fracture) or zero
!! (parametrized fracture).
!!
!! \note Only the 2D case is currently implemented.
!-----------------------------------------------------------------------
module mod_fracture_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_fracture_test_constructor, &
   mod_fracture_test_destructor,  &
   test_name, &
   test_description,&
   coeff_diff,&
   coeff_adv, &
   coeff_re,  &
   coeff_f,   &
   coeff_dir, &
   coeff_neu, &
   coeff_rob, &
   coeff_alpha

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 character(len=*), parameter ::    &
   test_name = "fracture"

! Module variables

 ! public members
 character(len=100), protected ::    &
   test_description(2)
 ! private members
 real(wp), parameter :: &
   ! permeability coefficients
   kf1 = 1.0_wp,    &
   kf2 = 2.0e-3_wp, &
   ! fracture dimension
   fd  = 0.01_wp,   &
   !fdg = fd,        & ! resolved fracture
   fdg = 0.0_wp,    & ! modeled fracture
   fh  = 0.8_wp
 integer :: ndim
 logical :: &
   mod_fracture_test_initialized = .false.
 character(len=*), parameter :: &
   this_mod_name = 'mod_fracture_test'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_fracture_test_constructor(d)
  integer, intent(in) :: d

  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_fracture_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   if(d.eq.2) then
     ndim = d
   else
     call error(this_sub_name,this_mod_name, &
                'This is a two-dimensional test case.')
   endif

   write(test_description(1),'(a)') &
     'Test case as in [Martin, Jaffre, Roberts, 2005]'
   write(test_description(2),'(a,e9.3,a,e9.3,a)') &
     '  permeabilities kf1 = ',kf1,', kf2 = ',kf2,'.'

   mod_fracture_test_initialized = .true.
 end subroutine mod_fracture_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_fracture_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_fracture_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_fracture_test_initialized = .false.
 end subroutine mod_fracture_test_destructor

!-----------------------------------------------------------------------
 
 pure function coeff_diff(x) result(mu)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: mu(size(x,1),size(x,1),size(x,2))

  integer :: l, i
  real(wp) :: xx, yy, kk

   mu = 0.0_wp
   do l=1,size(x,2)

     xx = x(1,l)
     yy = x(2,l)

     if(abs(xx).gt.fdg/2.0_wp) then ! outside fracture
       kk = 1.0_wp
     else ! inside fracture
       if(abs(yy).gt.fh/2.0_wp) then ! lateral regions
         kk = kf1
       else
         kk = kf2
       endif
     endif

     do i=1,size(x,1)
       mu(i,i,l) = kk
     enddo

   enddo
 
 end function coeff_diff
 
!-----------------------------------------------------------------------
 
 pure function coeff_adv(x) result(a)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: a(size(x,1),size(x,2))

   a = 0.0_wp
 end function coeff_adv
 
!-----------------------------------------------------------------------

 pure function coeff_re(x) result(sigma)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: sigma(size(x,2))
 
   sigma = 0.0_wp
 end function coeff_re
 
!-----------------------------------------------------------------------
 
 pure function coeff_f(x) result(f)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: f(size(x,2))

   f = 0.0_wp
 end function coeff_f
 
 
!-----------------------------------------------------------------------

 pure function coeff_dir(x,breg) result(d)
  real(wp), intent(in) :: x(:,:)
  integer, intent(in) :: breg
  real(wp) :: d(size(x,2))

   bregcase: select case(breg)
    case(2,8)
     d = 0.0_wp
    case(4,6)
     d = 1.0_wp
    case default
     d = huge(x(1,1))
   end select bregcase
 
 end function coeff_dir
 
!-----------------------------------------------------------------------

 pure function coeff_neu(x,breg) result(h)
  real(wp), intent(in) :: x(:,:)
  integer, intent(in) :: breg
  real(wp) :: h(size(x,2))

   bregcase: select case(breg)
    case(1,2,3,5,6,7)
     h = 0.0_wp
    case default
     h = huge(x(1,1))
   end select bregcase
 
 end function coeff_neu
 
!-----------------------------------------------------------------------

 pure function coeff_rob(x,breg) result(g)
  real(wp), intent(in) :: x(:,:)
  integer, intent(in) :: breg
  real(wp) :: g(size(x,2))

   g = 0.0_wp
 end function coeff_rob
 
!-----------------------------------------------------------------------

 pure function coeff_alpha(x) result(alpha)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: alpha(size(x,2))

  real(wp), parameter :: toll = 1.0e-5_wp
  integer :: l

   alpha = 0.0_wp
   do l=1,size(alpha)
     if( abs(x(1,l)).lt.toll ) then
       if( abs(x(2,l)).gt.fh/2.0_wp ) then
         alpha(l) = fd/(2.0_wp*kf1)
       else
         alpha(l) = fd/(2.0_wp*kf2)
       endif
     endif
   enddo

 end function coeff_alpha
 
!-----------------------------------------------------------------------

end module mod_fracture_test

